## Package to install: 
install.packages("remotes")
library(remotes)
remotes::install_version("DIDmultiplegtDYN",version="1.0.8",dependencies=TRUE)

# Libraries used
library(devtools)
options(rgl.useNULL=TRUE)
library(DIDmultiplegtDYN)
library(dplyr)
library(ggplot2)
library(plm)
library(xtable)
library(stringr)
library(fixest)

# Set working directory as appropriate 
setwd("/Users/michaelalbertus/Dropbox/Research/Indigenous Communities in Peru/Paper on Recognition and Conflict/JCR/R&R/Replication files")


############################################

## Table 1 prep

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

## set treatment and outcome
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# drop provinces without conflict/survey
df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)


######################################
## Table 1: Average Total Effects ####
  
effects=20
placebo=10
estt2.1<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant2.1 = mean(estt2.1$save_sample$outcome)
outsdt2.1 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt2.1$save_sample)))

effects=11
placebo=6
estt2.2<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)                        
outmeant2.2 = mean(estt2.2$save_sample$outcome)
outsdt2.2 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt2.2$save_sample)))


##########################
## Table 2: Robustness ###

## Model 1

# Drop 1980
df_analysis<-subset(df_analysis,df_analysis$p_year!=1980)

effects=11
placebo=5
estt3.1<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)  
outmeant3.1 = mean(estt3.1$save_sample$outcome)
outsdt3.1 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt3.1$save_sample)))
                                             

## Model 2

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)

df_analysis$largecomm = ifelse(df_analysis$families>2000,1,0)
df_analysis$largecomm = ifelse(is.na(df_analysis$largecomm),0,df_analysis$largecomm)
df_analysis<-subset(df_analysis,df_analysis$largecomm==0) 

effects=11
placebo=5
estt3.2<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant3.2 = mean(estt3.2$save_sample$outcome)
outsdt3.2 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt3.2$save_sample)))
               
## Model 3

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)

# Drop communities in the Amazon region
df_analysis<-subset(df_analysis,df_analysis$natural_regions!=2)

effects=11
placebo=6
estt3.3<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant3.3 = mean(estt3.3$save_sample$outcome)
outsdt3.3 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt3.3$save_sample)))
               

## Model 4

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$PROVINCIA,df_analysis$DISTRITO,sep="")
df_analysis$Comm_IDn = as.numeric(df_analysis$Comm_ID)

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(deptprovdist2) %>%
  dplyr::mutate(death_district = sum(comm_num_death_peace_full))%>%
  filter(death_district != 0)

effects=11
placebo=6
estt3.4<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant3.4 = mean(estt3.4$save_sample$outcome)
outsdt3.4 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt3.4$save_sample)))
               
                         
## Model 5

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$PROVINCIA,df_analysis$DISTRITO,sep="")
df_analysis$Comm_IDn = as.numeric(df_analysis$Comm_ID)

## set tx and outcome
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_mue

df_analysis<-df_analysis %>%
  dplyr::group_by(deptprovdist2) %>%
  dplyr::mutate(death_district = sum(comm_num_death_mue))%>%
  filter(death_district != 0)

effects=10
placebo=6
estt3.5<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant3.5 = mean(estt3.5$save_sample$outcome)
outsdt3.5 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt3.5$save_sample)))
               
## Model 6

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$PROVINCIA,df_analysis$DISTRITO,sep="")
df_analysis$Comm_IDn = as.numeric(df_analysis$Comm_ID)

## set tx and outcome
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_event_actos

df_analysis<-df_analysis %>%
  dplyr::group_by(deptprovdist2) %>%
  dplyr::mutate(death_district = sum(comm_num_event_actos))%>%
  filter(death_district != 0)

effects=10
placebo=6
estt3.6<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant3.6 = mean(estt3.6$save_sample$outcome)
outsdt3.6 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt3.6$save_sample)))
               
               
#####################################
## Table 3: Property Mechanism ######

## Model 1

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-subset(df_analysis,df_analysis$land_acquisition==3|df_analysis$land_acquisition==4)

effects=12
placebo=5
estt4.1<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant4.1 = mean(estt4.1$save_sample$outcome)
outsdt4.1 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt4.1$save_sample)))

## Model 2
 
######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-subset(df_analysis,df_analysis$land_acquisition==1)

effects=9
placebo=6
estt4.2<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant4.2 = mean(estt4.2$save_sample$outcome)
outsdt4.2 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt4.2$save_sample)))

## Model 3
 
######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# Create land pressure measure
df_analysis$pressure1 = df_analysis$pop_1961v2/df_analysis$tarea_ha

# subset to high pressure areas
df_analysis<-subset(df_analysis,!is.na(df_analysis$pressure1)&df_analysis$pressure1>median(df_analysis$pressure1,na.rm=T))

effects=10
placebo=7
estt4.3<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant4.3 = mean(estt4.3$save_sample$outcome)
outsdt4.3 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt4.3$save_sample)))

## Model 4
 
######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# Create land pressure measure
df_analysis$pressure1 = df_analysis$pop_1961v2/df_analysis$tarea_ha

# subset to low pressure areas
df_analysis<-subset(df_analysis,!is.na(df_analysis$pressure1)&df_analysis$pressure1<median(df_analysis$pressure1,na.rm=T))

effects=10
placebo=5
estt4.4<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant4.4 = mean(estt4.4$save_sample$outcome)
outsdt4.4 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt4.4$save_sample)))


## Model 5

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

## set tx and outcome
df_analysis$treatment<-df_analysis$treat_title
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)

effects=6
placebo=6
estt4.5<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_recog"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)                        
outmeant4.5 = mean(estt4.5$save_sample$outcome)
outsdt4.5 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_recog | Comm_ID+department^p_year,estt4.5$save_sample)))


######################################################
## Table 4: State-Community Relations Mechanism ######

## Model 1

maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$province,df_analysis$DISTRITO,sep="")
df_analysis$deptprovdist3 = toupper(df_analysis$deptprovdist2)

# Merge in district-level data
districtdata = read.csv("districtdata.csv")
df_analysis = merge(df_analysis,districtdata,by.x="deptprovdist3",by.y="deptprovdi2",all.x=T)

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
  
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# subset to places high road density
df_analysis<-subset(df_analysis,!is.na(df_analysis$d_roads_1973_all)&df_analysis$d_roads_1973_all>median(df_analysis$d_roads_1973_all,na.rm=T))

effects=11
placebo=6
estt5.1<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant5.1 = mean(estt5.1$save_sample$outcome)
outsdt5.1 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt5.1$save_sample)))


## Model 2

maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$province,df_analysis$DISTRITO,sep="")
df_analysis$deptprovdist3 = toupper(df_analysis$deptprovdist2)

# Merge in district-level data
districtdata = read.csv("districtdata.csv")
df_analysis = merge(df_analysis,districtdata,by.x="deptprovdist3",by.y="deptprovdi2",all.x=T)

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# subset to places with low road density
df_analysis<-subset(df_analysis,!is.na(df_analysis$d_roads_1973_all)&df_analysis$d_roads_1973_all<median(df_analysis$d_roads_1973_all,na.rm=T))

effects=11
placebo=6
estt5.2<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant5.2 = mean(estt5.2$save_sample$outcome)
outsdt5.2 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt5.2$save_sample)))


## Model 3
 
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$civilian_violence_peace_full

effects=11
placebo=6
estt5.3<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant5.3 = mean(estt5.3$save_sample$outcome)
outsdt5.3 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt5.3$save_sample)))
               
## Model 4

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$guerrilla_violence_peace_sl_mrta_full

effects=11
placebo=6
estt5.4<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)
outmeant5.4 = mean(estt5.4$save_sample$outcome)
outsdt5.4 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt5.4$save_sample)))


## Model 5

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$state_violence_peace_full

effects=11
placebo=6
estt5.5<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions","previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)


## Model 6

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

## set tx and outcome
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$guerrilla_violence_peace_sl_mrta_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
# Restrict to 1989-2000
df_analysis<-subset(df_analysis,df_analysis$p_year>1988)

effects=7
placebo=2
estt5.6<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)  
outmeant5.6 = mean(estt5.6$save_sample$outcome)
outsdt5.6 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt5.6$save_sample)))


## Model 7

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$state_violence_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
# Restrict to 1984-2000
df_analysis<-subset(df_analysis,df_analysis$p_year>1983)

effects=10
placebo=3
estt5.7<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)  
outmeant5.7 = mean(estt5.7$save_sample$outcome)
outsdt5.7 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt5.7$save_sample)))


######################################################
## Table 5: Participation in Rondas Campesinas #######

# Read in districts panel data 
df_district = read.csv("districtsrondas.csv")

## Model 1

estt6.1 <- feols(mronda_camp94 ~ sharerec_pcdistrict_pre1994post1980 + elev1k + slope + population_1972_10k + landreformarea + marxistleft + d_roads_1973_all_1k + shareCC| department,se="cluster",df_district)
summary(estt6.1)

## Model 2

estt6.2 <- feols(mronda_camp94 ~ sharerec_pcdistrict_pre1994post1980 + elev1k + slope + population_1972_10k + landreformarea + marxistleft + d_roads_1973_all_1k + shareCC +totalevents | department,se="cluster",df_district)
summary(estt6.2)

## Model 3

estt6.3 <- feols(mronda_camp94 ~ sharerec_pcdistrict_pre1994post1980 + elev1k + slope + population_1972_10k + landreformarea + marxistleft + d_roads_1973_all_1k +shareCC +totalevents| department,subset = df_district$population_1972<25000,se="cluster",df_district)
summary(estt6.3)

## Model 4

estt6.4 <- feols(mronda_camp94 ~ sharerec_pcdistrict_pre1994post1980 + elev1k + slope + population_1972_10k + landreformarea + marxistleft + d_roads_1973_all_1k + shareCC +totalevents| department,subset = df_district$sum_tarea/df_district$sarea_ha>.3,se="cluster",df_district)
summary(estt6.4)

## Model 5

estt6.5 <- feols(mronda_camp94 ~ sharerec_pcdistrict_pre1994post1980 + elev1k + slope + population_1972_10k + landreformarea + marxistleft + d_roads_1973_all_1k + shareCC +totalevents| department,subset = df_district$department!="Cajamarca",se="cluster",df_district)
summary(estt6.5)


########################
## Descriptive stats of CVR conflict events and deaths for footnote in paper

maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$PROVINCIA,df_analysis$DISTRITO,sep="")

df_analysis<-df_analysis %>%
  dplyr::group_by(deptprovdist2) %>%
  dplyr::mutate(death_district = sum(comm_num_event_actos))%>%
  filter(death_district != 0)

df_analysis$treatment<-df_analysis$treat_recog

df_analysis_d<-df_analysis %>%
   dplyr::group_by(Comm_ID) %>%
   dplyr::mutate(treat_comm_count = sum(treatment))%>%
   filter(treat_comm_count != 21)

summary(df_analysis_d$comm_num_event_actos)
sd(df_analysis_d$comm_num_event_actos)

summary(df_analysis_d$comm_num_death_mue)
sd(df_analysis_d$comm_num_death_mue)


###############################################################
####### Appendix


######################################
### Table VI.i. Clustering SEs at higher levels of aggregation 

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)


df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$PROVINCIA,df_analysis$DISTRITO,sep="")

effects=20
placebo=10
estt1.5<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T,cluster=c("deptprovdist2"))
outmeant1.5 = mean(estt1.5$save_sample$outcome)
outsdt1.5 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt1.5$save_sample)))


effects=11
placebo=6
estt1.6<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T,cluster=c("deptprovdist2"))                        
outmeant1.6 = mean(estt1.6$save_sample$outcome)
outsdt1.6 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estt1.6$save_sample)))



######################################
### Table VI.ii. History of prior community mobilization

## Model 1

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# subset to places with no recorded prior movements
df_analysis<-subset(df_analysis,is.na(df_analysis$movementcom)&is.na(df_analysis$movementdis)&is.na(df_analysis$movementprov)&is.na(df_analysis$petition1920s))

effects=11 
placebo=6 
estsamplesplit1.1<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)  
outmeant2.1 = mean(estsamplesplit1.1$save_sample$outcome)
outsdt2.1 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estsamplesplit1.1$save_sample)))
estsamplesplit1.1

## Model 2

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)
    
df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

# subset to places with prior movements
df_analysis<-subset(df_analysis,df_analysis$movementcom==1|df_analysis$movementdis==1|df_analysis$movementprov==1|df_analysis$petition1920s==1)

effects=13 
placebo=7 
estsamplesplit1.1<-did_multiplegt_dyn(df_analysis, outcome = "outcome", group = "Comm_ID", time = "p_year", 
                         treatment = "treatment", graph_off = TRUE, 
                         effects = effects, placebo = placebo, 
                         controls = c("has_open_concessions", "previous_mining_claims","treat_title"),
                         trends_nonparam = c("department"),ci_level=90,save_sample=T)  
outmeant2.1 = mean(estsamplesplit1.1$save_sample$outcome)
outsdt2.1 = sd(resid(feols(outcome ~ has_open_concessions+previous_mining_claims+treat_title | Comm_ID+department^p_year,estsamplesplit1.1$save_sample)))
estsamplesplit1.1


######################################
### Figures VI.iv. Event-study plots by perpetrator

## Civilian-perpetrated deaths
plt <- estt5.3$plot 
plt$layers[[1]] <- NULL 
plt$layers[[1]]$aes_params$colour <- c(rep("#00BFC4", 11), rep("#F7736B", 6)) 
plt$layers[[2]]$aes_params$colour <- c(rep("#00BFC4", 12), rep("#F7736B", 6)) 
plt <- plt + geom_hline(yintercept = 0, linetype = "dashed", color = "black") +  
  theme(plot.title = element_text(hjust = 0.5), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black")) + 
  scale_x_continuous(breaks=seq(-10,20,5)) +  
  ylab("Estimates") + ggtitle("") 
plt

## Guerrilla-perpetrated deaths
plt <- estt5.4$plot 
plt$layers[[1]] <- NULL 
plt$layers[[1]]$aes_params$colour <- c(rep("#00BFC4", 11), rep("#F7736B", 6)) 
plt$layers[[2]]$aes_params$colour <- c(rep("#00BFC4", 12), rep("#F7736B", 6)) 
plt <- plt + geom_hline(yintercept = 0, linetype = "dashed", color = "black") +  
  theme(plot.title = element_text(hjust = 0.5), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black")) + 
  scale_x_continuous(breaks=seq(-10,20,5)) +  
  ylab("Estimates") + ggtitle("") 
plt

## State-perpetrated deaths
plt <- estt5.5$plot 
plt$layers[[1]] <- NULL 
plt$layers[[1]]$aes_params$colour <- c(rep("#00BFC4", 11), rep("#F7736B", 6)) 
plt$layers[[2]]$aes_params$colour <- c(rep("#00BFC4", 12), rep("#F7736B", 6)) 
plt <- plt + geom_hline(yintercept = 0, linetype = "dashed", color = "black") +  
  theme(plot.title = element_text(hjust = 0.5), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black")) + 
  scale_x_continuous(breaks=seq(-10,20,5)) +  
  ylab("Estimates") + ggtitle("") 
plt


#####################
## Section V. Descriptive stats of where recognition occurs as a function of easy/hard places

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$deptprovdist2 = paste(df_analysis$DEPARTAMENTO,df_analysis$province,df_analysis$DISTRITO,sep="")
df_analysis$deptprovdist3 = toupper(df_analysis$deptprovdist2)

districtdata = read.csv("districtdata.csv")
df_analysis = merge(df_analysis,districtdata,by.x="deptprovdist3",by.y="deptprovdi2",all.x=T,no.dups=T)
df_analysis$priormovement<-ifelse(df_analysis$movementcom==1|df_analysis$movementdis==1|df_analysis$movementprov==1|df_analysis$petition1920s==1,1,0)
df_analysis$priormovement<-ifelse(is.na(df_analysis$priormovement),0,df_analysis$priormovement)

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)

## Drop to sample of CCs that change recognition status and do cross-tabs ####
df_analysis_d<-df_analysis %>%
   dplyr::group_by(Comm_ID) %>%
   dplyr::mutate(treat_comm_count = sum(treatment))%>%
   filter(treat_comm_count != 21)

df_analysis_mayors = subset(df_analysis_d,df_analysis_d$p_year==1993)
df_analysis_priormovement = subset(df_analysis_d,df_analysis_d$p_year==2000)

table(df_analysis_mayors$norenewalvacanciapost,df_analysis_mayors$treatment)
table(df_analysis_priormovement$priormovement,df_analysis_priormovement$treatment)

# For Appendix: t-test difference in recognition rates in places with mayors vs without 
df_analysis_mayors$recby2000 = ifelse((!is.na(df_analysis_mayors$FECHA_year)&df_analysis_mayors$FECHA_year<2001)|(!is.na(df_analysis_mayors$FECHA.RECONOCIMIENTO_year)&df_analysis_mayors$FECHA.RECONOCIMIENTO_year<2001),1,0)
table(df_analysis_mayors$norenewalvacanciapost,df_analysis_mayors$recby2000)
t.test(df_analysis_mayors$recby2000 ~ df_analysis_mayors$norenewalvacanciapost,data=df_analysis_mayors)

#######################
## Section V: Line plots of trends in recognition and violence

## First in early 1980s Emergency Zones
comm_conflict<-read.csv("CCpanel.csv")

# This is the EZ 1984 group at the cradle of the conflict. See p. 357 of this: http://lanic.utexas.edu/project/laoap/desco/desco00002.pdf
comm_conflict = subset(comm_conflict,(comm_conflict$department=="AYACUCHO"&province!="PARINACOCHAS"&province!="PAUCAR")|(comm_conflict$department=="HUANCAVELICA")|(comm_conflict$department=="APURIMAC"&(comm_conflict$province=="ANDAHUAYLAS"|comm_conflict$province=="CHINCHEROS")))
comm_conflict = subset(comm_conflict,comm_conflict$p_year<1996)

#line plot: sum number of community deaths per year 
comm_conflict_deaths_by_year <- comm_conflict %>% group_by(p_year) %>% summarize(sum = sum(comm_num_death_peace_full)) 

comm_confl_deaths <- ggplot(comm_conflict_deaths_by_year, aes(x=p_year, y=sum, group =1))+
  geom_line(color = "darkred", linewidth = .6)+
  labs(x="Year", y="Conflict Deaths in Communities")+theme(aspect.ratio = 1/2)

ggsave("comm_confl_deaths_EZ1984.pdf", width=1600, height=800, units="px")

#sum number of recognitions per year, cumulative 
comm_recog_by_year <- comm_conflict %>% group_by(p_year) %>% summarize(sum = sum(treat_recog)) 

comm_recogs <- ggplot(comm_recog_by_year, aes(x=p_year, y=sum, group =1))+
  geom_line(color = "darkred", linewidth = .6)+
  labs(x="Year", y="Recognized Communities")+theme(aspect.ratio = 1/2)

ggsave("comm_recogs_EZ1984.pdf", width=1600, height=800, units="px")

## Next examine timing in Ayacucho, Huancavelica, Junin
comm_conflict<-read.csv("CCpanel.csv")

## Select these one at a time and rerun code on each to generate plots
comm_conflict<-subset(comm_conflict,comm_conflict$department=="AYACUCHO")
#comm_conflict<-subset(comm_conflict,comm_conflict$department=="HUANCAVELICA")
#comm_conflict<-subset(comm_conflict,comm_conflict$department=="JUNIN")

#line plot: sum number of community deaths per year 
comm_conflict_deaths_by_year <- comm_conflict %>% group_by(p_year) %>% summarize(sum = sum(comm_num_death_peace_full)) 

comm_confl_deaths <- ggplot(comm_conflict_deaths_by_year, aes(x=p_year, y=sum, group =1))+
  geom_line(color = "darkred", linewidth = .6)+
  labs(x="Year", y="Conflict Deaths in Communities")+theme(aspect.ratio = 1/2)

ggsave("comm_confl_deaths_dept.pdf", width=1600, height=800, units="px")

#sum number of recognitions per year, cumulative 
comm_recog_by_year <- comm_conflict %>% group_by(p_year) %>% summarize(sum = sum(treat_recog)) 

comm_recogs <- ggplot(comm_recog_by_year, aes(x=p_year, y=sum, group =1))+
  geom_line(color = "darkred", linewidth = .6)+
  labs(x="Year", y="Recognized Communities")+theme(aspect.ratio = 1/2)

ggsave("comm_recogs_dept.pdf", width=1600, height=800, units="px")


######################
# Section VII: Descriptive stats

######## ---------------- read in data ---------------- ########
maindata<-read.csv("CCpanel.csv")

df_analysis<-maindata

df_analysis$treatment<-df_analysis$treat_recog
df_analysis$outcome<-df_analysis$comm_num_death_peace_full

df_analysis<-df_analysis %>%
  dplyr::group_by(province) %>%
  dplyr::mutate(death_province = sum(comm_num_death_peace_full))%>%
  filter(death_province != 0)

df_analysis_d<-df_analysis %>%
   dplyr::group_by(Comm_ID) %>%
   dplyr::mutate(treat_comm_count = sum(treatment))%>%
   filter(treat_comm_count != 21)

summary(df_analysis_d$treatment)
sd(df_analysis_d$treatment)
summary(df_analysis_d$outcome)
sd(df_analysis_d$outcome)
summary(df_analysis_d$has_open_concessions)
sd(df_analysis_d$has_open_concessions)
summary(df_analysis_d$previous_mining_claims)
sd(df_analysis_d$previous_mining_claims)
summary(df_analysis_d$treat_title)
sd(df_analysis_d$treat_title)

df_analysis_d$largecomm = ifelse(df_analysis_d$families>2000,1,0)
df_analysis_d$largecomm = ifelse(is.na(df_analysis_d$largecomm),0,df_analysis_d$largecomm)
summary(df_analysis_d$largecomm)
sd(df_analysis_d$largecomm)
df_analysis_d$amazon = ifelse(df_analysis_d$natural_regions==2,1,0)
summary(df_analysis_d$amazon)
sd(df_analysis_d$amazon,na.rm=T)

df_analysis_d$insecurepr<-ifelse(df_analysis_d$land_acquisition==3|df_analysis_d$land_acquisition==4,1,0)
summary(df_analysis_d$insecurepr)
sd(df_analysis_d$insecurepr,na.rm=T)

summary(df_analysis_d$civilian_violence_peace_full)
sd(df_analysis_d$civilian_violence_peace_full)
summary(df_analysis_d$guerrilla_violence_peace_sl_mrta_full)
sd(df_analysis_d$guerrilla_violence_peace_sl_mrta_full)
summary(df_analysis_d$state_violence_peace_full)
sd(df_analysis_d$state_violence_peace_full)

# Land pressure measure
df_analysis_d$pressure1 = df_analysis_d$pop_1961v2/df_analysis_d$tarea_ha
summary(df_analysis_d$pressure1)
sd(df_analysis_d$pressure1,na.rm=T)

# State presence via roads
df_analysis_d$deptprovdist2 = paste(df_analysis_d$DEPARTAMENTO,df_analysis_d$province,df_analysis_d$DISTRITO,sep="")
df_analysis_d$deptprovdist3 = toupper(df_analysis_d$deptprovdist2)
# Merge in roads data
districtdata = read.csv("districtdata.csv")
df_analysis_d = merge(df_analysis_d,districtdata,by.x="deptprovdist3",by.y="deptprovdi2",all.x=T)
roadstats = summary(df_analysis_d$d_roads_1973_all,digits=6)
statsdf = as.data.frame(as.list(roadstats))
statsdf
sd(df_analysis_d$d_roads_1973_all,na.rm=T)


